La plupart des difficultés d’installation des librairies venaient de uwot, dont une version précédente était nécessaire et le couple SeuratWrappers et monocle3 dont leur installation demandent un token. Afin d’obtenir le token, je vous conseille de regarder ce guide https://gist.github.com/z3tt/3dab3535007acf108391649766409421 . Une fois le token bien enregistré dans .Renviron et R redémarré, la série suivante d’installations devrait permettre l’installation des packages problématiqes:
install.packages("devtools")
library(devtools)
devtools::install_version("uwot", version = "0.1.10", repos = "http://cran.us.r-project.org")
library(uwot)
install.packages("Seurat")
library(Seurat)
BiocManager::install("Rsamtools")
install.packages("R.utils")
library(Rsamtools)
library(R.utils)
remotes::install_github('satijalab/seurat-wrappers')
library(SeuratWrappers)
devtools::install_github('cole-trapnell-lab/leidenbase')
BiocManager::install("limma")
BiocManager::install("batchelor")
library(limma)
librairy(batchelor)
devtools::install_github('cole-trapnell-lab/monocle3')
library(monocle3)
Pour toutes les libraries restantes, un simple install.packages("...") suffira.
#devtools::install_version("uwot", version = "0.1.10", repos = "http://cran.us.r-project.org")
library(uwot)
library(Seurat)
library(tximport)
library(ggplot2)
library(corrplot)
library(network)
library(Signac)
library(SeuratWrappers)
library(monocle3)
library(Matrix)
library(patchwork)
library(gridExtra)
set.seed(1234)
setwd("~/mydatalocal/scarabi/results")
Récupération des matrices des 7 échantillons:
sampsf <- c("SRR8257100","SRR8257101","SRR8257102","SRR8257103","SRR8257104","SRR8257105","SRR8257106")
files <- file.path(
paste("~/mydatalocal/scarabi/data/quant/",sampsf,"/alevin/quants_mat.gz", sep=""))
txis <- lapply(files, function(f) tximport(files = f, type="alevin"))
Création des objects Seurat associés :
seu_objs <- lapply(seq_along(txis), function(i){
# min.cells = 3, min.features = 200 from Ryu et al.
s <- CreateSeuratObject(counts = txis[[i]]$counts , min.cells = 3, min.features = 200, project = sampsf[i])
})
Pour l’instant, on le récupère que les samples du Wild Type (non mutants):
scarabWT <- merge(x = seu_objs[[1]], y = unlist(seu_objs[2:4], recursive = F), add.cell.ids = sampsf[1:4])
Ensuite, on étudie nos échantillons afin de ne conserver qu’une partie de nos données, en enlevant les données ne nous semblant pas significatives:
scarabWT[["percent.mt"]] <- PercentageFeatureSet(scarabWT, pattern = "ATM")
scarabWT[["percent.chloro"]] <- PercentageFeatureSet(scarabWT, pattern = "ATC")
quant <- quantile(scarabWT[["percent.mt"]]$percent.m,0.95 )
scarabWT <-subset(scarabWT, subset = percent.mt < quant & percent.chloro < 0.2)
Voici le résultat obtenu :
VlnPlot(scarabWT, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.chloro"), ncol = 4)
plot1 <- FeatureScatter(scarabWT, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(scarabWT, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))
Ensuite, on normalise nos données et l’on garde qu’un certain nombre nombre (pour ma version, 2000) des cellules les plus intéressantes:
scarabWT <- NormalizeData(scarabWT, normalization.method = "LogNormalize", scale.factor = 10000)
nb_cells = 2000
scarabWT <- FindVariableFeatures(scarabWT, selection.method = "vst", nfeatures = nb_cells)
Maintenant, on lance une ACP sur nos données :
all.genes <- rownames(scarabWT)
scarabWT <- ScaleData(scarabWT, features = all.genes)
scarabWT <- RunPCA(scarabWT, features = VariableFeatures(object = scarabWT))
Et l’on obtient ceci :
VizDimLoadings(scarabWT, dims = 1:2, reduction = "pca")
DimPlot(scarabWT, reduction = "pca")
DimHeatmap(scarabWT, dims = 1:15, cells = 500, balanced = TRUE)
L’étude des composantes principales a pour but de déterminer la dimension de notre problème. Afin de pouvoir déterminer le plus de types cellulaires possibles dans nos données, il est préférable d’avoir une dimension élevée. Cependant, chaque nouvelle composante principale doit aussi avoir un intérêt à discriminer nos données et doit donc être significative. Alors, pour se faire, on utilise deux analyses supplémentaires : JackStraw et l’ElbowPlot. Voici les graphiques obtenues :
scarabWT <- JackStraw(scarabWT, num.replicate = 100)
scarabWT <- ScoreJackStraw(scarabWT, dims = 1:20)
JackStrawPlot(scarabWT, dims = 1:15, xmax = 1, ymax = 1)
ElbowPlot(scarabWT)
Ce que l’on y obverse sur l’ElbowPlot est à partir de la 10ème composante principale, les suivantes ne sont plus si intéressantes que ça. On a donc décidé (arbitrairement) de se restreindre à une dimension de 10.
Afin de générer la carte UMAP (avec RunUMAP), on a besoin au préalable de convertir nos données en clusters, via un graphe pondéré en calculant les voisins de chaque cellules.
scarabWT <- FindNeighbors(scarabWT, dims = 1:10)
scarabWT <- FindClusters(scarabWT, resolution = 0.8)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 11321
## Number of edges: 342018
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9284
## Number of communities: 21
## Elapsed time: 1 seconds
scarabWT <- RunUMAP(scarabWT, dims = 1:5, return.model = TRUE, umap.method = "uwot-learn")
On obtient alors les résultats suivants :
p8 <-DimPlot(scarabWT, reduction = "umap")
p8
ggsave(p8,file = "~/mydatalocal/scarabi/results/image/UMAP.png",width=30,height=15,units="cm")
Afin d’annoter la carte UMAP en associant les différentes régions avec des types cellulaires, nous avons utilisé deux méthodes.
Le but de cette méthode est d’utiliser des marqueurs génétiques pour nos différents types cellulaires.
Markers = read.csv("../data/Markers.csv", sep = '\t', h=T)
Table = table(Markers$Preferential.expression.in.root)
Voici le tableau représant le nombre de marqueurs par type cellulaire :
Table
##
## Cortex
## 12
## Endodermis
## 9
## Epidermis - Non-hair Cells (Atrichoblasts)
## 11
## Epidermis - Root Hair Cells (Trichoblasts)
## 10
## Epidermis and Lateral Root Cap
## 5
## Meristem (Young) Endodermis/Cortex
## 7
## Root Cap (Columella)
## 17
## Stele (Young Vascular Tissue)
## 13
Ensuite, pour chacun des types cellulaires, on affiche sur la carte UMAP seulement l’expression des gènes de ce type cellulaire.
Markers$Locus<-gsub(" ","",Markers$Locus)
Markers$Preferential.expression.in.root<-gsub("/"," ",Markers$Preferential.expression.in.root)
lm<-split(Markers,Markers$Preferential.expression.in.root)
system("mkdir -p ~/mydatalocal/scarabi/results/image")
system("rm -r ~/mydatalocal/scarabi/results/image/*")
output <- lapply(names(lm),function(x){f<-FeaturePlot(scarabWT, features = lm[[x]]$Locus)
ggsave(f,file=paste0("~/mydatalocal/scarabi/results/image/",x,".png"),width=40,height=40,units="cm")
})
datascore <- data.frame(lapply(names(lm),function(x){score=colMeans(scarabWT@assays$RNA[lm[[x]]$Locus,])
}))
names(datascore)<-make.names(names(lm))
scarabWT <- AddMetaData(scarabWT, metadata = datascore)
g<-FeaturePlot(scarabWT, features=names(datascore))
g
ggsave(g,file = "~/mydatalocal/scarabi/results/image/type_cellulaire.png",width=40,height=40,units="cm")
Via des prélévements on a aussi accès à l’expression des gènes par type cellulaire
samps<-read.table("~/mydatalocal/scarabi/data/Counts_Salmon/metadata_Li2016.txt", sep="\t")
#On trie les samples par types cellulaires
samps <- samps[order(samps$V3),]
#On supprime les types cellulaires non déterminants
samps <- samps[!samps$V3%in%c("whole root","whole root 1","cycloheximide treatment","cycloheximide mock"),]
ech<-samps$V1
Ensuite, pour chacun de ces samples on récupère les matrices de compte et via à une table, on réassocie ensemble les noms des gènes.
files <- file.path(
paste("~/mydatalocal/scarabi/data/Counts_Salmon/",ech,"/quant.sf", sep=""))
files<-files[file.exists(files)]
tx2gene<-read.table("~/mydatalocal/scarabi/processed_data/txp2gene.tsv")
names(tx2gene)<-c("TXNAME","GENEID")
tx2gene<-unique(tx2gene)
head(tx2gene)
## TXNAME GENEID
## 1 AT1G01010.1 AT1G01010
## 7 AT1G01020_P2 AT1G01020
## 15 AT1G01020_P6 AT1G01020
## 21 AT1G01020_P1 AT1G01020
## 30 AT1G01020_P3 AT1G01020
## 38 AT1G01020_P4 AT1G01020
On peut alors récupérer les données et on obtient alors la matrice de compte (normalisée) suivante :
txis <- lapply(files, function(f) {
tab<- tximport(files = f, type="salmon", tx2gene=tx2gene)
return(tab$abundance)
})
tabpur<-as.data.frame(txis)
ech2=sapply(files,function(f){strsplit(f,"/")[[1]][6]})
#Changer le nom des colonnes du tableau
names(tabpur)<-make.names(ech2)
head(tabpur)
## SRR3664382 SRR3664383 SRR3664385 SRR3664395 SRR3664406 SRR3664417
## AT1G01010 1.716811 1.896926 1.545711 15.843074 20.268393 27.714170
## AT1G01020 7.369509 11.111502 9.280125 17.836873 19.710903 14.096914
## AT1G01030 0.987394 1.474313 0.275222 3.689367 3.250175 3.115996
## AT1G01040 9.021629 8.665490 4.201520 5.220210 4.976830 6.818153
## AT1G01046 0.000000 0.000000 1.405640 0.000000 0.000000 0.000000
## AT1G01050 80.577198 96.215644 71.923622 80.756158 70.374419 67.141841
## SRR3664396 SRR3664397 SRR3664398 SRR3664419 SRR3664420 SRR3664421
## AT1G01010 1.871991 3.618997 0.644235 11.358003 13.906525 8.793016
## AT1G01020 18.363256 16.313037 7.771199 9.363583 9.615792 10.181512
## AT1G01030 0.227622 0.047495 0.050757 0.094970 0.213155 0.238876
## AT1G01040 6.728216 5.175364 2.719236 5.278344 5.429668 5.547827
## AT1G01046 0.000000 0.000000 0.630508 0.000000 0.000000 0.819966
## AT1G01050 57.092115 71.099493 33.008388 92.716870 85.178594 103.502638
## SRR3664402 SRR3664403 SRR3664404 SRR3664428 SRR3664435 SRR3664436
## AT1G01010 23.222836 5.580218 9.812438 3.259710 3.046452 4.292104
## AT1G01020 5.655989 10.831351 7.492238 3.077802 2.955736 1.872432
## AT1G01030 6.632473 3.191752 3.568764 0.050170 0.000000 0.234531
## AT1G01040 8.307703 7.156862 9.470113 2.188932 1.946090 1.110394
## AT1G01046 0.000000 0.000000 1.257230 0.000000 0.000000 0.000000
## AT1G01050 28.472910 46.433552 36.656115 39.044528 34.581970 41.300804
## SRR3664405 SRR3664407 SRR3664408 SRR3664422 SRR3664423 SRR3664424
## AT1G01010 1.595829 1.669632 1.277080 30.262334 25.156696 24.216769
## AT1G01020 9.652126 6.374632 8.140137 7.261154 4.492351 6.229815
## AT1G01030 0.271967 0.308613 0.746134 0.150571 0.042560 0.035130
## AT1G01040 3.506512 3.004830 3.717348 5.564117 5.734301 5.665490
## AT1G01046 0.000000 0.548677 0.000000 0.000000 0.000000 0.000000
## AT1G01050 59.719122 50.988552 59.998486 48.434264 35.376794 36.978165
## SRR3664374 SRR3664375 SRR3664437 SRR3664376 SRR3664377 SRR3664378
## AT1G01010 15.959187 18.168734 16.162023 57.228042 93.080597 83.976778
## AT1G01020 3.445357 3.176128 2.993439 8.075620 10.339301 7.353513
## AT1G01030 1.131589 1.439347 1.989967 3.381933 7.935872 5.471184
## AT1G01040 2.952577 2.444447 2.646486 9.455829 11.725951 11.487360
## AT1G01046 0.000000 0.000000 0.000000 1.121207 1.812935 6.681024
## AT1G01050 40.327258 25.524147 34.209081 30.594107 28.945098 25.546687
## SRR3664389 SRR3664391 SRR3664425 SRR3664426 SRR3664427 SRR3664379
## AT1G01010 4.049538 9.052426 2.310489 1.541885 2.024963 4.263242
## AT1G01020 6.610719 3.517385 17.869868 16.397251 21.811314 5.658542
## AT1G01030 0.176140 0.000000 0.291288 0.322293 0.495996 0.358554
## AT1G01040 18.197588 7.864732 3.984504 3.071081 3.582128 2.855099
## AT1G01046 15.835439 0.614266 0.349925 0.000000 0.000000 0.000000
## AT1G01050 75.779156 59.369302 132.357691 101.977525 91.140975 57.979350
## SRR3664380 SRR3664381 SRR3664372 SRR3664373 SRR3664384 SRR3664386
## AT1G01010 3.593588 5.582260 64.629978 9.211267 39.910340 55.078668
## AT1G01020 6.424684 5.219054 6.489138 2.905927 13.995652 6.056696
## AT1G01030 0.267133 0.488625 0.188868 0.658552 0.420462 0.000000
## AT1G01040 2.350481 2.728850 14.268800 6.399822 16.890305 15.307501
## AT1G01046 0.000000 0.000000 0.584967 0.000000 2.268211 0.000000
## AT1G01050 51.868707 48.501560 48.210757 35.352159 73.961593 43.024911
## SRR3664387 SRR3664388 SRR3664392 SRR3664393 SRR3664394 SRR3664412
## AT1G01010 53.791810 37.618442 9.763332 18.649549 5.724260 0.604879
## AT1G01020 1.283277 5.794505 19.191938 8.081400 8.419938 17.998218
## AT1G01030 0.074136 0.000000 0.223327 0.069081 0.177317 0.622511
## AT1G01040 20.716409 17.283298 6.199051 5.724517 7.265130 4.821658
## AT1G01046 0.000000 0.721670 0.000000 0.748399 0.000000 0.000000
## AT1G01050 55.027919 43.526574 58.147394 42.592574 49.841970 48.721253
## SRR3664413 SRR3664414
## AT1G01010 2.017441 2.845328
## AT1G01020 14.953029 16.290637
## AT1G01030 0.362785 1.603667
## AT1G01040 4.778282 8.002214
## AT1G01046 0.000000 0.000000
## AT1G01050 43.536514 50.942496
Enfin, pour chaque cluster, on regarde l’expression moyenne des gènes dans le cluster. Puis on ne garde que les gènes en commun, on calcule la matrice de corrélation entre les clusters et les samples et finalement on change le nom des samples (aux colonnes) par le type cellulaire associé. On obtient alors ceci :
avg.e <- AverageExpression(scarabWT)
scarabWT_avg=data.frame(avg.e)
genes_scarabi <- rownames(scarabWT_avg)
genes_li <- rownames(tabpur)
genes_common <- genes_scarabi[genes_scarabi%in%genes_li]
countsLi_norm_avg_alt_sum_c <- tabpur[genes_common,]
scarabWT_avg_c <- scarabWT_avg[genes_common,]
corLi_scarab_spearman <- cor(scarabWT_avg_c,countsLi_norm_avg_alt_sum_c,method="spearman")
colnames(corLi_scarab_spearman) <- lapply(colnames(corLi_scarab_spearman), function(name){samps[samps$V1==name,3]})
corrplot(corLi_scarab_spearman, method="color", is.corr=F, tl.col = as.color(colnames(corLi_scarab_spearman)) )
Alors, pour chaque cluster, en prennant par exemple le type cellulaire avec la plus grande expression, on peut annoter la carte UMAP :
cluster.id <- max.col(corLi_scarab_spearman)
cluster.id <- colnames(corLi_scarab_spearman)[cluster.id]
names(cluster.id) <- levels(scarabWT)
scarabWT <- RenameIdents(scarabWT, cluster.id)
p7 <- DimPlot(scarabWT, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
p7
ggsave(p7,file = "~/mydatalocal/scarabi/results/image/UMAP_label.png",width=15,height=15,units="cm")
Afin de comparaitre nos mutants avec le Wild Type, on recommence toutes les opérations jusqu’à réobtenir la carte UMAP :
files <- file.path( paste("~/mydatalocal/scarabi/data/quant/SRR8257105/alevin/quants_mat.gz", sep=""))
txi <- tximport(files = files, type="alevin")
scarabmut <- CreateSeuratObject(counts = txi$counts , min.cells = 3, min.features = 200, project = sampsf[6])
scarabmut[["percent.mt"]] <- PercentageFeatureSet(scarabmut, pattern = "ATM")
scarabmut[["percent.chloro"]] <- PercentageFeatureSet(scarabmut, pattern = "ATC")
quant <- quantile(scarabmut[["percent.mt"]]$percent.m,0.95 )
scarabmut <-subset(scarabmut, subset = percent.mt < quant & percent.chloro < 0.2)
scarabmut <- NormalizeData(scarabmut, normalization.method = "LogNormalize", scale.factor = 10000)
scarabmut <- FindVariableFeatures(scarabmut, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(scarabmut)
scarabmut <- ScaleData(scarabmut, features = all.genes)
scarabmut <- RunPCA(scarabmut, features = VariableFeatures(object = scarabmut))
scarabmut <- FindNeighbors(scarabmut, dims = 1:10)
scarabmut <- FindClusters(scarabmut, resolution = 0.8)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2371
## Number of edges: 67404
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8897
## Number of communities: 14
## Elapsed time: 0 seconds
scarabmut <- RunUMAP(scarabmut, dims = 1:5, return.model = TRUE, umap.method = "uwot-learn")
p6 <-DimPlot(scarabmut, reduction = "umap")
p6
ggsave(p6,file = "~/mydatalocal/scarabi/results/image/UMAP_mut.png",width=30,height=15,units="cm")
avg.e <- AverageExpression(scarabmut)
scarabmut_avg=data.frame(avg.e)
genes_scarabi <- rownames(scarabmut_avg)
genes_li <- rownames(tabpur)
genes_common <- genes_scarabi[genes_scarabi%in%genes_li]
countsLi_norm_avg_alt_sum_c <- tabpur[genes_common,]
scarabmut_avg_c <- scarabmut_avg[genes_common,]
corLi_scarab_spearman <- cor(scarabmut_avg_c,countsLi_norm_avg_alt_sum_c,method="spearman")
colnames(corLi_scarab_spearman) <- lapply(colnames(corLi_scarab_spearman), function(name){samps[samps$V1==name,3]})
corrplot(corLi_scarab_spearman, method="color", is.corr=F, tl.col = as.color(colnames(corLi_scarab_spearman)) )
cluster.id2 <- max.col(corLi_scarab_spearman)
cluster.id2 <- colnames(corLi_scarab_spearman)[cluster.id2]
names(cluster.id2) <- levels(scarabmut)
scarabmut <- RenameIdents(scarabmut, cluster.id2)
p5 <- DimPlot(scarabmut, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
p5
ggsave(p5,file = "~/mydatalocal/scarabi/results/image/UMAP_mut_Label.png",width=30,height=15,units="cm")
Maintenant que l’on a les cartes UMAP des deux, un moyen intéressant de les comparer est de projeter les cellules et cluster du mutant sur la carte du Wild Type. On obtient alors ceci :
scarab.anchors <- FindTransferAnchors(reference = scarabWT, query = scarabmut,
dims = 1:30, reference.reduction = "pca")
scarabmut <- MapQuery(anchorset = scarab.anchors, reference = scarabWT, query = scarabmut,
reference.reduction = "pca", reduction.model = "umap")
p1 <- DimPlot(scarabWT, reduction = "umap", group.by = "seurat_clusters", label = TRUE, label.size = 3,
repel = TRUE) + NoLegend() + ggtitle("Reference annotations")
p2 <- DimPlot(scarabmut, reduction = "ref.umap", group.by = "seurat_clusters", label = TRUE,
label.size = 3, repel = TRUE, split.by = "orig.ident") + NoLegend() + ggtitle("Query transferred labels")
p1 + p2
ggsave(p1+p2,file = "~/mydatalocal/scarabi/results/image/projection_mutant.png",width=30,height=15,units="cm")
Comme on peut l’observer, toute la partie poilue a disparu chez le mutant. Alors, une analyse supplémentaire est possible : on peut étudier le processus de différentiation des cellules de ce groupe via le PseudoTime.
atricho.cds <- scarabWT[, scarabWT$seurat_clusters %in% c(15,11,6,9,17)]
atricho.cds <- as.cell_data_set(atricho.cds)
atricho.cds <- cluster_cells(cds = atricho.cds, reduction_method = "UMAP")
atricho.cds <- learn_graph(atricho.cds, use_partition = TRUE)
##
|
| | 0%
|
|======================================================================| 100%
# meristem to root
meristemcells <- names(scarabWT$seurat_clusters[scarabWT$seurat_clusters==17])
atricho.cds <- order_cells(atricho.cds, reduction_method = "UMAP", root_cells = meristemcells[22:22])
scarabWT <- AddMetaData(
object = scarabWT,
metadata = atricho.cds@principal_graph_aux@listData$UMAP$pseudotime,
col.name = "Trichoblast"
)
p=FeaturePlot(scarabWT, c("Trichoblast"), pt.size = 0.1)
p3 <- plot_cells(
cds = atricho.cds,
color_cells_by = "pseudotime",
label_branch_points=FALSE,
show_trajectory_graph = TRUE,
)
p3
grid.arrange(p, p3, ncol=2)
ggsave(grid.arrange(p, p3, ncol=2),file="~/mydatalocal/scarabi/results/image/peudoTime.png",width = 10, height = 5)